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Abstract 

The main objective of this work is to explore aspects of stochastic resonance (SR) in noisy 
bistable, symmetric systems driven by subthreshold periodic rectangular external signals possessing 
a large duty cycle of unity. Using a precise numerical solution of the Langevin equation, we carry 
out a detailed analysis of the behavior of the first two cumulant averages, the correlation function 
and its coherent and incoherent parts. We also depict the non-monotonic behavior versus the 
noise strength of several SR quantifiers such as the average output amplitude, i.e. the spectral 
amplification (SPA), the signal-to-noise ratio (SNR) and the SR-gain. In particular, we find that 
with subthreshold amplitudes and for an appropriate duration of the pulses of the driving force 
the phenomenon of stochastic resonance (SR), is accompanied by SR-gains exceeding unity. This 
analysis thus sheds new light onto the interplay between nonlinearity and the nonlinear response 
which in turn yields nontrivial, unexpected SR-gains above unity. 

PACS numbers: 05.40.-a, 05.10.Gg, 02.50.-r 



' ticasado@us.es| |http : //numerix . us . es| 



1 



I. INTRODUCTION 



Most of the studies on the phenomenon of Stochastic Resonance (SR) in dynamical 



ml m aynam 



for 



systems have been devoted to systems driven by sinusoidal terms (see 
reviews). Several analytical approximations have been put forward to explain SR. In the 
approach of McNamara and Wiesenfeld the Langevin dynamics is replaced by a reduced 
two-state model that neglects the intra-well dynamics. The general ideas of Linear Response 
Theory (LRT) have been applied to situations where the input amplitude is small 



1 neory [LtiL ^ 
HQ- In 0, 



111 ] the Floquet theory has been applied to the corresponding Fokker-Planck 



description. For very low input frequencies, an adiabatic ansatz has been invoked . Even 
though these alternative analytical approaches provide an explanation of SR for different 
regions of parameter values, their precise limits of validity remain to be determined. In 



recent work, Il2 



13l | we have explored the validity of LRT for sinusoidal and multifrequency 



input signals with low frequency. Our results indicate a breakdown of the LRT description 
of the average behavior for low frequency, subthreshold amplitude inputs. 

Several quantifiers have been used to characterize SR in noisy, continuous, systems. The 
average output amplitude, or the spectral amplification (SPA), has been studied in Refs. ^, 



111 ] and the phase of the output average in Refs. [Ij,ll5|, respectively. Those parameters as 



well as the signal-to-noise ratio (SNR) |6|, exhibit a non-monotonic behavior with the noise 
strength which is representative of SR. An important quantity is the SR-gain, defined as 
the ratio of the SNR of the output over the input SNR. It has been repeatedly pointed out 
that the SR-gain can not exceed unity as long as the system operates in a regime described 



by LRT 
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j. Beyond LRT there exists no physical reason that prevents the SR-gain to 
be larger than 1, as it has been demonstrated in for super-threshold sinusoidal input 
signals, and in analog experiments in 
Fourier components and a small duty cycle 



2G| for subthreshold input signals with many 

12|. 

In this work, we will make use of numerical solutions of the Langevin equation following 
the methodology in fl^ , to analyze SR in noisy bistable systems, driven by periodic piecewise 
constant signal with two amplitude values of opposite signs (rectangular signal) (see Fig. HI). 
There are several relevant time scales in the dynamics of these systems: i) tasym, the time 
interval within each half period of the driving force, during which the diffusing particle 
sees an asymmetric constant two-well potential; ii) tinter, the time scale associated with the 
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inter- well transitions in both directions; and iii) tmtra, the time scale associated to intra- well 
dynamics. The inter-well and intra-well time scales depend basically on the noise strength 
D and the amplitude of the driving term. The dependence of these two time scales with 
those parameters is certainly very different, being more pronounced for titter- Typically, for 
the range of parameter values associated with SR, the intra-well time scale is shorter than 
the inter-well one. 

We will evaluate the long-time average behavior of the output and the second cumulant. 



These two quantities were studied some time ago by two of us [2l| for periodic rectangular 
driving signals. Here, we will further extend our work to the analysis of the correlation 
function and its coherent and incoherent parts. The knowledge of all these quantities pro- 
vides a very useful information for the explanation of SR as indicated by the non-monotonic 
behavior with the noise strength of the output amplitude and the SNR. In particular, the 
knowledge of the incoherent part of the correlation function is of outmost importance for a 
correct determination of the SNR. Furthermore, for a given subthreshold amplitude, we will 
demonstrate that, if there exist a range of noise values such that tenter is shorter than tasym, 
then it is possible to observe stochastic amplification and, simultaneously, SR-gains larger 
;han unity. This is strictly forbidden by linear response theory, as we have recently shown 



17[. Thus, a simultaneous appearance of stochastic amplification and SR-gains above 1 
implies a strong violation of linear response theory. 



II. MODEL SYSTEM AND SR QUANTIFIERS 

Let us consider a system characterized by a single degree of freedom, x, subject to the 
action of a zero average Gaussian white noise with {C,{t)C,{s)) = 2D6{t — s) and driven by 
an external periodic signal F{t) with period T. In the Langevin description, its dynamics is 
generated by the equation 

x{t) = -U'[xit)]+Fit)+at). (1) 

The corresponding linear Fokker-Planck equation (FPE) for the probability density P{x, t) 
reads 

%P{x,t) = t{t)P{x,t), (2) 



3 



where 



U'{x)-F{t) + D^ 



(3) 



dx 

In the expressions above, U'{x) represents the derivative of the potential U{x). In this 
work, we will consider a bistable potential U{x) = — + The periodicity of the 

external driving F{t) allows its Fourier series expansion in the harmonics of the fundamental 
frequency Q = 2tt/T, i.e., 

oo 

F{t) = [fn cos(nfit) + gn sin(nfit)] , (4) 

n=l 

with the Fourier coefficients, and Qn, given by 

fn = ^ r dtF{t)cosinnt), 
I Jo 

gn = ^ f dtF{t)sm{nQt). (5) 
I Jo 

Here, we are assuming that the cycle average of the external driving over its period equals 
zero. In this work, we will focus our attention to multi-frequency input forces with "rectan- 
gular" shape given by 

f A : < t < ^ 
[-A : ^<t<T 

as sketched in Fig. ^ The external force remains constant at a value A during each half- 
period and changes sign for the second half of the period. The duty cycle of the signal is 
defined by the time span the signal is nonzero over the total period of the signal; thus, the 
rectangular signal in Eq. (jH)) consequently possesses a duty cycle of unity. 

The two-time correlation function {x{t + r)x(t))oo in the limit t oo is given by 

/oo /"OO 
dx' x'Pooix',t) / dxxPi\i{x,t + T\x',t), (7) 
-oo J —oo 

where P^{x,t) is the time-periodic, asymptotic long time solution of the FPE and the 
quantity Pi|i(x, t + r|x',t) denotes the two-time conditional probability density that the 
stochastic variable will have a value near x at time t + r if its value at time t was exactly 
x' . It can been shown ^ that, in the limit t cxo, the two-time correlation function 
{x{t + r)x(t))oo becomes a periodic function of t with the period of the external driving. 
Then, we define the one-time correlation function, C(r), as the average of the two-time 
correlation function over a period of the external driving, i.e., 

C{r) = - [ dt(x(t + r)x(t))oo. (8) 
I Jo 



The correlation function C(r) can be written exactly as the sum of two contributions: a 



coherent part, Ccoh{T), which is periodic in r with period T, and an inco 



lerent part which 



decays to for large values of r. The coherent part Ccohi^) is given by 

Ccohir) = ^ r At {x{t + r))oo(a;(t))oo, (9) 



TJo 

where (x(t))oo is the average value evaluated with the asymptotic form of the probability 
density, Pqo{x, t). 

According to McNamara and Wiesenfeld 0|, the output SNR is defined in terms of the 
Fourier transform of the coherent and incoherent parts of C(r). As the correlation function 
is even in time and we evaluate its time dependence for r > 0, it is convenient to use its 
Fourier cosine transform, defined as 

~ 2 /""^ ~ 

C{uo) = - AtC{t)cos{ujt)- C{t)= AujC{uo)cos{uot). (10) 

TT JO JO 

The value of the output SNR is then obtained from: 

^^^^ ^ lim,^o.^"+;d^(7H _ ^^^^ 

C incoh (^) 



mingfn 

m 



rom the same contri- 
The periodicity of 



Note that this definition of the SNR differs by a factor 2, stemrnin^ 
bution at = — fi, from the definitions used in earlier works 
the coherent part gives rise to delta peaks in the spectrum. Thus, the only contribution to 
the numerator in Eq. (fTT|) stems from the coherent part of the correlation function. The 
evaluation of the SNR requires the knowledge of the Fourier components of Ccohi'T) and 
Cincohij) at the fundamental frequency of the driving force. Thus, rather than knowledge of 
the entire Fourier spectrum, only two well defined numerical quadratures are needed. These 
are: 

Rout = 77^' (12) 



where 



2 f'^ 

Q« = ^ dr Ccoh(r) cos(i7r), (13) 



and 

2 /'°° 

Qi = - drCincohir) cosinr). (14) 

TT JO 

The SNR for an input signal F{t) + ^{t) is given by 



The SR-gain G is consequently defined as the ratio of the SNR of the output over the SNR 
of the input, namely, 

G = (16) 



III. NUMERICAL SOLUTION 

Stochastic trajectories, x^^\t), are generated by numerically integrating the Langevin 
equation [Eq. (^] for every realization j of the white noise ^(t), starting from a given initial 
condition xn. The numerical solution is based on the algorithm developed by Greenside and 
Helfand (22, |23| (consult also the Appendix in Ref. 17[). After allowing for a relaxation 
transient stage, we start recording the time evolution of each random trajectory for many 
different trajectories. Then, we construct the long time average value, 

N 



1 



and the second cumulant, 

1 ^ 

{x\t))^-{xit))l = -Y.[x^'\t)_ 



(17) 



N 



(18) 



where is the number of stochastic trajectories considered. We also evaluate the two-time 
{t and r) correlation function, i.e.. 



N 



(x(t + r)x(t))c 
as well as the product of the averages 
{x{t + r))oo(a;(t))oo = 



(19) 



1 



N 



1 



N 



— Yx^^^t) 



(20) 



The correlation function C(r) and its coherent part Ccohi'^) are then obtained using their 
definitions in Eqs. (jHl) and Q, performing the cycle average over one period of t. The differ- 
ence between the values of C(r) and CcohiT) allows us to obtain the values for Cincohij)- It 
is then straightforward to evaluate the Fourier component of Ccohi'^) and the Fourier trans- 
form of Cincohij) at the driving frequency by numerical quadrature. With that information, 
the numerator and the denominator for the output SNR [cf. Eqs. (O, (fTSj) and (HH)], as 
well as the SR-gain [cf. Eq. |TB|]. are obtained. 



IV. RESULTS 



A. Response to a rectangular driving force with fundamental frequency 17 = 0.01 

Consider an external driving of the type sketched in Fig. with parameter values 
Q = 0.01, A = 0.25. This amplitude is well below its threshold value defined, for each 
driving frequency, as the minimum amplitude that can induce repeated transitions between 
the minima of U{x) in the absence of noise. For the input considered here, the threshold 
amplitude is At ~ 0.37. Note that this threshold value for the amplitude increases with 
increasing driving frequency. 

In Fig. 121 we depict with several panels the behavior of the first two cumulants, {x{t))oo 
and (x^(t))oo — {x(t))l^, for several representative values of D [from top to bottom D = 0.02 
(panel a), D = 0.04 (panel h), D = 0.06 (panel c), D = 0.1 (panel d) and D = 0.2 (panel 
e)]. Notice that due to the transients, the time at which we start recording data, t = in 
the graphs, does not necessarily coincide with the start of an external cycle. The average is 
periodic with the period of the driving force, while the second cumulant, due to the reflection 
symmetry of the potential il2], is periodic with a period half of the period of the forcing 
term. 

Next we consider the case of small noise intensity D (say, D = 0.02 as in panel a). The 
noise induces jumps between the wells. In each random trajectory, a jump between the 
wells has a very short duration, but the instants of time at which they take place for the 
different stochastic trajectories are randomly distributed during a half-cycle. At this small 
noise strength, the jumps are basically towards the lowest minimum. Thus, because of this 
statistical effect, the average behavior depicts the smooth evolution depicted in panel a of 
Fig. 121 without sudden transitions between the wells. The evolution of the second cumulant 
adds relevant information. The fact that it is rather large during most of a period indicates 
that the probability density Poo{x,t) is basically bimodal during most of the external cycle. 
It is only during very short time intervals around each half-period that the probability 
distribution becomes monomodal around one of the minima, and, consequently, the second 
cumulant is small. The bimodal character arises from the fact that the noise is so small in 
comparison with the barrier heights that jumps over the barrier are rather infrequent during 
each half a period. 
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As the noise strength increases, the time evolution of (x(t))oo follows closely the shape 
of the external force, (cf. see panels b, c, d in Fig. |21for 0.04 < D < 0.1). This behavior 
indicates that, for these parameter values, the jumps in the different random trajectories are 
concentrated within short time intervals around the instants of time at which the driving 
force switches sign. The second cumulant remains very small during most of a period, except 
for short time intervals around the switching times of the external driver. Thus, for these 
intermediate values of D, the probability distribution, Poo{x,t) is basically monomodal, 
except for small time intervals around the switching instants of time of the periodic driver. 
Finally, as the noise strength is further increased {D > 0.1), the probability distribution 
remains very broad most of the time. Even though a large majority of random trajectories 
will still jump over the barrier in synchrony with the switching times of the driver, the 
noise is so large in comparison with the barrier heights, that the probability of crossings 
over the barrier in both directions can not be neglected at any time during each half cycle. 
The probability distribution remains bimodal during a whole period, but asymmetric: the 
larger fraction of the probability accumulates around the corresponding lower minima of the 
potential during each cycle. Therefore, the average output amplitude decreases, while the 
second cumulant depicts plateaus at higher values than for smaller noise strengths (compare 
panels e and c in Fig. |21). 

In Fig. El we plot the coherent (left panels), Cco/i(t), and incoherent (right panels), 
Cincoh{T)i components of the correlation function C(r) for the same parameter values as in 
Fig. 121 The coherent part shows oscillations with a period equal to that of the driving force. 
Its shape changes with D. The amplitude of the coherent part does not grow monotonically 
with D. Rather, it maximizes at ~ 0.06, which is consistent with the observed behavior 
of {x(t))oo in Fig. |21 This is expected as the evaluation of Cco/i(t) involves only the time 
behavior of {x{t))oQ at two different instants of time. 

Two features of the behavior of Cmcohi'T) are relevant: its initial value and its decay time. 
The initial value of the incoherent contribution, Cj„co/i(0), is given by the cycle average of 
the second cumulant. It has a non-monotonic behavior with D. For D = 0.02, Cj„coh(0) is 
large, consistent with the fact that the second cumulant at this noise strength is appreciably 
different from during a substantial part of a period. As the value of D increases, {D < 0.1), 
the value of Cincohi^) decreases. This is expected as the second cumulant is large just during 
those small time intervals where most of the forward transitions take place every half-period. 
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For still larger values of D there are frequent forward and backward jumps that keep the 
stationary probability bimodal, and therefore, the initial value Cmcohi^) increases. For 
D = 0.02, the decay of Cincoh{T) is very slow, although the decay time is still shorter than 
the duration of half a period of the driving force. As D increases, the decay time of the 
incoherent part becomes shorter. It is worth to point out that the intra-well noisy dynamics 
manifests itself in the behavior of Cincohi'^)- This is most clearly confirmed noticing the fast 
initial decay observed in panel d. For smaller values of D, this feature is masked by the long 
total relaxation time scale, while for D = 0.2, the noise strength is so large that there is not 
a clear-cut separation between inter and intra-well time scales. 

The above considerations allow us to rationalize the behavior of the several quantifiers 
used to characterize SR. Their behaviors with D for A = 0.25, Q = 0.01 are depicted in 
Fig. m It should be noticed that the lowest value of the noise strength used in the numerical 
solution of the Langevin equation is -D = 0.01. For this noise strength, the values of Qu 
and Rout are very small, although not zero. For even lower noise strengths the task becomes 
computationally very demanding and expensive, due to the extremely slow decay of the 
correlations. For D sufficiently small, however, one does expect Qu to be larger than Qi, 
and, consequently, an increase of the numerical Rout as D is lowered. 

The quantity Qu defined in Eq. i\Vd\i depicts a non-monotonic dependence on D typical 
of the SR phenomenon. Its behavior is expected from the dependence of the amplitude of 
Ccohir) with D in Fig. H 

A non-monotonic behavior with D for the numerically evaluated Qi is also observed. The 
initial value, Cincoh{0) and the decay time of CmcohiT) are important in the evaluation of 
Qi [see Eq. fll4|) ]. For D = 0.01, the decay time of the incoherent part of the correlation 
function is longer than half a period of the driving force, while for D = 0.02, it is somewhat 
shorter than T/2. Consistently with Eq. ()14j) . the value of the integral for D = 0.01 is 
smaller than for D = 0.02. As D is further increased, the infiuence of the cosine factor in 
Eq. (fT^ becomes less important as the decay time is much shorter than T/2. The drastic 
fall in the Qi values observed for 0.02 < D < 0.1 is due to the decrease of Cincoh{0) with D 
(see panels a-c in Fig. ^ and the shortening of the decay time. As D is increased further, 
Cincohi^) increases and, consequently, Qi also increases slightly. 

Taking into account the definition of the SNR, [cf. Eq. ()12|) ]. its behavior with D is not 
surprising. The numerically obtained SNR peaks a.t D = 0.08, a slightly different value of 
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D from the one at which Qu peaks. 

B. Anomalous SR-gain behavior for subthreshold driving 

The SR-gain is defined in Eq. ()16|) . The numerically determined SR-gain shows a most 
interesting feature: we observe a non-monotonic behavior versus D, with values for the gain 
exceeding unity (!) for a whole range of noise strengths. This is strictly forbidden within 
LRT jl^; therefore, the fact that the SR-gain can assume values larger than unity refiects a 
manifestation of the inadequacy of LRT to describe the system dynamics for the parameter 
values considered. 

To rationalize this anomalous SR-gain behavior, we notice that the role of the noise in the 
dynamics is twofold. On the one hand, it controls the decay time of Cmcohi'^)- On the other 
hand, the noise value is relevant to ellucidate whether the one-time probability distribution 
is basically monomodal or bimodal during most of the cycle and, consequently, it controls 
the initial value Cmco?t(0). As discussed above, if D is small, the decay time is very large 
compared to tasym, and the one-time probability distribution is essentially bimodal. For large 
values of D, the decay of Cincoh{T) is fast enough, and the distribution is also bimodal. The 
large SR-gain obtained here requires the existence of a range of intermediate noise values 
such that: i) Cj„coh(T) decays in a much shorter time scale than tasym and, ii) the one-time 
probability distribution remains monomodal during most of the external cycle. 

C. Response to a rectangular input driver with fundamental frequency 0, = 0.1 

As mentioned before, there are several time scales which are important for the phe- 
nomenon of stochastic amplification and gain. In the previous subsection we have considered 
an input frequency small enough so that the inequality tasym > tmter holds for a range of 
noise values. Next, we shall analyze the system response to a driving force with fundamental 
frequency Q = 0.1, ten times larger than in the previous case. We will take the same input 
amplitude as in the previous case, A = 0.25, which is still subthreshold. For this input fre- 
quency, the threshold value for the amplitude strength is determined by numerically solving 
the deterministic equation to yield At ~ 0.39. 

The behavior of the first two cumulants for several values of the noise strength is depicted 
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in Fig. |31 [from top to bottom D = 0.02 (panel a), D = 0.04 (panel h), D = 0.06 (panel c), 
D = 0.1 (panel d) and D = 0.2 (panel e)]. For all values of D, the second cumulant remains 
large for most of each half-period. By contrast with the lower frequency case, we detect no 
values of D for which the probability distribution is monomodal for a significant fraction of 
each half period. 

In Fig. inithe behavior of the coherent (left panels) and incoherent (right panels) parts of 
the correlation function are presented for (from top to bottom) D = 0.02 (panel a), D = 0.04 
(panel h), D = 0.06 (panel c), D = 0.1 (panel d) and D = 0.2 (panel e). The amplitude 
of the coherent oscillations shows a nonmonotonic behavior with D. The incoherent part 
has initial values which remain very large in comparison with the corresponding ones for 
Q = 0.01 (compare with Fig. consistently with the large value of the second cumulant. 
The decay times are roughly the same for both frequencies. 

In Fig. [7|we show the behavior of the several SR quantifiers as a function of D. The 
comparison of Figs. IH and indicates that Qu , Qi and Rout have the same qualitative 
behavior for both frequencies. The non-monotonic dependence on D of Qu and Rout are 
indicative of the existence of SR (for both frequencies) for the subthreshold input amplitude 
and in the ranges of D values considered. The most relevant quantitative difference is that 
for = 0.1 the SR-gain remains less than unity. 



V. CONCLUSIONS 



With this work, we have analyzed the phenomenon of SR within the context of a noisy, 
bistable symmetric system driven by time periodic, rectangular forcing possessing a duty 
cycle of unity. The numerical solution of the Langevin equation allows us to analyze the 
long-time behavior of the average, the second cumulant and the coherent and incoherent 
parts of the correlation function. For subthreshold input signals we determined the SNR, 
together with its numerator and denominator evaluated separately, for a wide range of noise 
strengths D. 

As a main result we find the simultaneous existence of a typical non-monotonic behavior 
versus the noise strength D of several quantifiers associated to SR; in particular SR-gains 
larger than unity are possible for a subthreshold rectangular forcing possessing a duty cycle 
of unity. This finding is at variance with the recent claim in Refs. |2o[ 24] that pulse- 
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like signals with small duty cycles are needed to obtain SR-gains larger than unity. This 
unexpected result occurs indeed whenever the inequality tasym > Unter holds. This is most 
easily achieved with low frequency inputs. As the input frequency increases, that inequality 
is not satisfied for sufficiently small values of D: even though SR then still exists, it is not 
accompanied by SR-gains exceeding unity. Furthermore, in Refs. |2o[ 24 1, SR-gains larger 



than unity are only obtained with input amplitudes larger than 0.8 Ay. By contrast, in this 
work, we have shown that such a large value for the input amplitude is not needed (we have 
used A ^ 0.68 At). 

The simultaneous occurrance of SR and SR-gains larger than unity is associated to the 
fact that, for some range of noise values, the decay time of the incoherent part of the corre- 
lation function is much shorter than tasym and also the probability distribution is basically 
monomodal during most of the cycle of the driving force. 
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FIG. 1: Sketch of a rectangular periodic signal with duty cycle 1, amplitude A and period T 
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FIG. 2: Time behavior of the average {x{t))oo (solid lines) and the second cumulant {x^{t))oo — 
{x{t))1^ (dashed lines) for a rectangular driving force with duty cycle 1, fundamental frequency 
= 0.01 and subthreshold amplitude A = 0.25 for several values of the noise strength: D = 0.02 
(panel a), D = 0.04 (panel b), D = 0.06 (panel c), D = 0.1 (panel d), D = 0.2 (panel e). Notice 
that due to the transients, t = in the graphs, does not necessarily coincide with the start of an 
external cycle. 
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FIG. 3: Time behavior of CcohiT) (left panels) and Cincohi.^) (right panels) for a rectangular driving 
force with duty cycle 1, fundamental frequency 17 = 0.01 and subthreshold amplitude A = 0.25 for 
several values of the noise strength: D = 0.02 (panel a), D = 0.04 (panel h), D = 0.06 (panel c), 
D = 0.1 (panel d), D = 0.2 (panel e). 





FIG. 4: Dependence with D of several SR quantifiers: the numerator of the SNR {Qu), its 
denominator (Q/), the output SNR {Rout) and the SR-gain (G) for a rectangular driving force with 
duty cycle 1, fundamental frequency S7 = 0.01 and subthreshold amplitude A = 0.25. 
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FIG. 5: Same as in Fig. El but for n = 0.1 
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FIG. 6: Same as in Fig. 01 but for n = 0.1. 
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FIG. 7: Same as in Fig. H but for O = 0.1. 
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